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By implementing the exact density matrix for the rotating anisotropic harmonic trap, we derive 
a class of very fast and accurate fourth order algorithms for evolving the Gross-Pitaevskii equation 
in imaginary time. Such fourth order algorithms are possible only with the use of forward, positive 
time step factorization schemes. These fourth order algorithms converge at time-step sizes an 
order-of-magnitude larger than conventional second order algorithms. Our use of time-dependent 
factorization schemes provides a systematic way of devising algorithms for solving this type of 
nonlinear equations. 



I. INTRODUCTION 



The dynamics of a fast rotating Bose-Einstein condensate (BEC) has been studied extensively in terms of the Gross- 
Pitaevskii (GP) equation^. By evolving the GP equation in imaginary time, it is easy to determine the ground state 
properties of the condensate, such as the formation of vortex-arrays and giant vortices^. It has been known for 
some time that the first order pseudo-spectral, split-operator method^ is a very fast way of solving the non-linear 
Schrodinger equation. However, first or second order split operator (SO) methods^ and Crank-Nickolson (CN) 
algorithms withi or without splitting ignore the time-dependence of the non-linear potential and converge linearly or 
quadratically only at very small time steps. Bandaruk and Shen^ have applied higher order decomposition schemes 
with negative coefficients to solve the real time non-linear Schrodinger equation. Due to the difficulty of estimating 
the non-linear potential at intermediate time, they have not demonstrated that their higher order algorithms actually 
converge with accuracy beyond second order. In any case, their negative coefficient algorithms cannot be used for 
imaginary time evolution because negative time steps will result in an unbounded diffusion kerneU2iiiii2iI£. 

In this work, we derive a class of very accurate fourth order factorization algorithm for solving the GP equation in 
imaginary time. These algorithms are made possible by the confluence of three key ideas: 1) The density matrix for a 
rotating anisotropic harmonic oscillator can be solved exactly. 2) The time-dependence of the non-linear potential can 
be systematically accounted for in factorization algorithms. 3) Forward, all positive time step algorithms^i^iSiiLiSiiS 
are now available for solving imaginary time evolution equations. 

In the next section, we show how the density matrix of the harmonic oscillator can be exactly implemented as an 
algorithm. This obviates the need to expand in harmonic eigenstates^S. In Section III, by exact diagonalization, we 
generalize the result to the case of a rotating anisotropic harmonic trap. In Section IV, we describe the time-dependent 
form of the fourth-order forward algorithm for solving the GP equation. In Section V, we compare the convergence 
of various algorithms. We summarize our conclusions in Section VI. 



II. EXACT ALGORITHM FOR THE HARMONIC OSCILLATOR 



Consider the 1-D harmonic oscillator Hamiltonian operator given by 

H = T + V = ^p 2 + \jj 2 x 2 . (2.1) 
Its imaginary time propagator (or density matrix) can be exactly decomposed as 

e ~r{T+V) = e - TCv V e -rC T T e -r Cv V^ (2 2) 

where Cy and Ct are functions of r to be determined. To show this, we simply compare the matrix elements on both 
sides. For the RHS, we have (ignoring normalization factors) 

(x'\e- rCvV e- TC ^ T e- TCvV \x) = ^c v ^^ -,f % -,c v ^,\ (2 3) 
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For the LHS of (|2.2|l . the exact density matrix element is known^i 
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If we keep only the first term, we have a second order algorithm. Keeping the first two terms gives a fourth order 
algorithm, keeping the first three terms gives a sixth order algorithm, etc.. 

The exact factorization (|2.2I) is possible for the harmonic oscillator because its Hamiltonian is quadratic and higher 
order commutators are either zero or simply proportional to T or V. The harmonic oscillator is characterized by two 
key commutators, 



[V,[T,V]} = 2tu 2 V, 



(2- 



[T, [V,T]} = 2uj 2 T. 



(2.9) 



Because of these two equalities, all higher order commutators can be subsumed back to the original operators T and 
V. To see how this exact decomposition comes about, let's begin with the simple second order decomposition, 



e _i T y e _ T T e -I.y 



exp 
exp 



-r(T + V) + ^r 3 ([V, [T, V}} - 2[T, [V,T]}) + 0(r 5 ) 
-t{T + V) + ^t 3 (2oj 2 V - WT) + 0(r 5 ) 



(2.10) 



Since the error terms are proportional to the original operators, they can be symmetrically moved back to the LHS 
to yield, 



'•(l-a«V)V p -T(l+i«V)r.-T^-^T")V _ p -r(T+y)+0(r°) 



(2.11) 



The decomposition of the LHS is then correct to fourth order. The coefficients agree with the expansion (I2.6|) and 
(|2.7|) . This example makes it clear that the exact expansion only depends on the abstract commutator relations 12. 8|) 
and H2.9|) . and is independent of the specific representation of the 1-D harmonic oscillator. Also, if we exchange the 
operators T <-> V, the coefficients are unchanged. Thus we can also factorize exactly via 



e -r(T+V) = e ~rCvT ^-tCtV e ~rC v T _ 



For real time propagation, we only need to set r = it to get the corresponding coefficients, 

1 - cos(wi) sinfwi) 
Cv — — ; — : — - — — and Ct = 



(2.12) 



(2.13) 



tut sin(wt) ' tut 

For either real or imaginary time evolution, one iterates the discretized wave function forward in time via 

H>(t + At)} = e - AT(T+v) | V(r)) (2.14) 

If the exact density matrix l|2.4|l were used directly in coordinate space, that would incur a slow, N x N matrix multi- 
plication of the Gaussian kernel e 2tC t . The advantage of the factorized form is that this matrix multiplication 

can be avoided by going to k-space via FFT and multiplying the k-space wave function point-by-point by e _rC ' T 2' £ . 
This is then an order 7Vln2 N operation, much faster than the N x TV coordinate space matrix multiplication. 
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III. EXACT ALGORITHM FOR A ROTATING ANISOTROPIC HARMONIC TRAP 

Consider now the case of an rotating anisotropic harmonic potential with Hamiltonian 

H = \{pl+p 2 v ) + \u 2 x x 2 + ~u>y -Q{xp y - yp x ) . (3.1) 

This is a well-studied problem in nuclear physics^ 2 ". Its diagonalization is greatly simplified 2 ^ if we characterize the 
anisotropy via the deformation parameter S 

u 2 x = (1 + 6)uj 2 , u, y = (l-6)u 2 , (3.2) 

measure lengths in units of the oscillator length I = and express H and Q in units of loq. The resulting 

dimensionless Hamiltonian is then 

H = \ipl+P 2 y ) + ^(1 + S)x 2 + \{1 - S)y 2 - fi (xp y - yp x ) , (3.3) 
where O = £1/ujq. To diagonalize this Hamiltonian, we introduce two new sets of canonical variables, 

Qi = ai(cx - sp y ), Pi = — {cp x + sy), (3.4) 

Oil 

Q 2 = a 2 (cy - sp x ), P 2 = —(cp y + sx), (3.5) 

OL 2 

where oti are normalization constants, and c = cos(</>), s — sin(</>). One can check that the canonical commutator 
relations are indeed satisfied, 

[Q i ,P j ] = i5 ij . (3.6) 

In terms of {Qi, Pi}, because of the way we have parameterized the anisotropy and expressed everything in terms of 
ujo, the coefficients of both P2Q1 and P1Q2 can be made to vanish with a single condition: 

on 

tan(20) = — . (3.7) 


Using on to normalize the P 2 terms with unit coefficient, the resulting Hamiltonian can be written as 

H = T 1 + V 1 +T 2 + V 2 = X -Pl + Ul\Q\ + X -Pl + \vl\Q\, (3.8) 



where 



o 



r 2 = 1- 1 + 1^2 + 4^2 ) 



a 2 - 2 = l + i-^p + AW, (3.9) 



with 



Also, from 1)3. 7fl . we have 



Vl\ = 1 + fl 2 + y/S 2 + 4fi 2 , 

nl = 1 + n 2 - Vs 2 + m 2 . (3.10) 



2 s 2 = 1 
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£ 

2c 2 = 1 + - = . (3.11) 

Vs 2 + m 2 
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At 11 = 0, the phase angle <fi = 0. As increases, the phase angle approaches 45° asymptotically. Thus s and c 
in (|3.11l) are both positive. However, as increases, fi^ crosses zero and becomes negative at fi = \fl — 8. At this 
critical rotation rate, the Coriolis force overcomes the weaker harmonic potential in the y-direction and the anisotropic 
harmonic oscillator is unstable. 0§ emerges positive again when a^ 2 crosses zero and turn negative at f2 2 = 1 + 6. 
Thus f2| is negative over the interval 1 — 6 < fl 2 < 1 + S. This is an instability of the rotating harmonic oscillator, 
not necessary that of the Gross-Pitaevskii equation. We will come back to this point in Section VI. Note also that for 
5 = 0, the algorithm is stable up to CI = 1. 

Eq. 1|3.8[) consists of two independent harmonic oscillator with different frequency. The two exact algorithms must 
be applied in sequence. However, since T\ and V2 only depend on p x and y, they should be placed next to each other 
so that both can be evaluated in the same mixed representation described below. Similarly, T2 and V\ only depend 
on x and p y . We therefore use the following factorization for each algorithm, 

e _ r ( Tl +Vi) =e -rC v (l)T le -TC T (l)Vi e -TC v (l)ri ) ( 3 _ 12 ) 
e -T(T 2 +V 2 ) =e -rC v {2)V 2e -TC T (2)T 2e -TG v (2)V 2 ^ 

and interlaced them as follow: 

e -r(T 1 +V 1 +T 2 + V 2 ) _ e -TCv(l)T 1 ^rCv(2)V2 e -TCT(l)V 1 -TCT(2)T 2e -rCv(l)T 1 ^rCv(2)V2 (3 ^\ 

Here we use the shorthand notations CV(1) = Cy(J7i), CV(2) = CV(^2), etc.. To implement (|3.14f) . let 

4>(x, „) = -L / d Px 4,(p x , y) jv* x , (3.15) 



/2tt 

and 



i;(p x ,y) = -L, I </rr;r. ,/;■(■ \ (3.16) 



4>{x,p y ) = -= I dyip(x,y)e 



/2tt 
1 

2^ 



dyd Px ^(p x ,y)e ip * x - ip yy. (3.17) 



The operators T± and V2 are diagonal in the representation ip(Px, y) and T2 and V± are diagonal in the representation 
ip(x,p y ). In practice, tp(x,y) is discretized as an N x N complex array and its Fourier transform is computed using 
the discretized FFT. Thus the exact algorithm consists of four steps: 

1. Compute the forward N-1D transform ip{p x ,y) from ip{x,y) and multiply ^{p x ,y) grid-point by grid-point by 
e -TCV(i)Ti-rCV(2)V2 ; wnere T\ and V2 are now understood to be functions of p x and y. 

2. Compute the 2D transform i/j(x,p y ) from the updated ijj(Px,y) and multiply ij)(x,py) by e - TC T(i)Vi-rC T {2)T 2 ^ 
where V% and T2 are now functions of x and p y . 

3. Compute the inverse 2D transform from the updated ip{x,p v ) back to ip{Px,y) and multiply rp(p x ,y) by 

e -TCv(l)T 1 -rCv(2)V 2 _ 

4. Compute the backward N-1D transform from the updated ip(p x ,y) back to ^(x,y). 

Thus the algorithm can be implemented with only three 2D-FFT. (One 2D-transform = 2N ID-transforms.) This is 
only one 2D-FFT more than solving the non-rotational case. 



IV. SOLVING THE GROSS-PITAEVSKII EQUATION 



Denoting now the entire rotating trap Hamiltonian l|3.8|l as the operator 

T = Ti + Vx+T 2 + V2, 



(4.1) 
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the corresponding 2D Gross-Pitaevskii equation is 

(T + g\iP\ 2 )iP(x,y) =^(x,y). (4.2) 
The condensate ground state can be projected out by imaginary time evolution: 

ip oc lim V(r) = lim e~ T [ T + y ( r )] +T ^/;(0). (4.3) 

r — >oc r — >oo 

The chemical potential // is determined by preserving the wave function's normalization to unity. This will be taken 
for granted and this term will be ignore in the following discussion. Since tp(r) is time-dependent, we have explicitly 
indicated that the Gross-Pitaevskii potential 

V{t)= 9 \^{t)\\ (4.4) 

is also time-dependent. 

In general, to solve (|4.3ll by factorization algorithms, one must apply rules of time-dependent factorizationiS* 2 ^: 

the time- dependent potential must be evaluated at an intermediate time equal to the sum of time steps of all the T 
operators to its right. For example, the first order algorithm 1A is 

V^Ar) = e- ATT e- ATy (°^(0) (4.5) 

and the first order algorithm IB is 

V-(Ar) = e - ATV(A ^ e- ATT V>(0). (4.6) 

While algorithm 1A is straightforward, IB requires that the potential be determined from the wave function to be 
computed. This self-consistency condition can be solved by iterative methods described below. 

In contrast to real time propagation, the wave function in imaginary time converges quickly to an approximate 
ground state depending on At and produces a ip(Ar) that differs from tp(0) only by a normalization constant. Thus 
after some initial iterations, the normalized g\ip(AT)\ 2 is independent of r and can be replaced by g\ip(0)\ 2 . This 
replacement can be justified only at small At when the approximate wave function is close to the exact ground state 
and is unchanging in time. (For real time propagation, the wave function is always changing with time, and one 
cannot justify this replacement even at small At.) At larger At, the approximate ground state may not be a discrete 
bound state and the algorithm may fail catastrophically. Thus if one approximates g\tp(AT)\ 2 by c?|i/)(0) | 2 in l|4.fi[l . 
then the algorithm is still first order, but only at very small At. We shall refer to this version of the algorithm as 
1B0. 

Wc define the second order algorithm 2A as 

^(At) = e-* A ^ AT )e~ ATT e-^ A ^(°V(0) (4.7) 

and algorithm 2B as 

V>(At) = e-5 ATT e- ATy ( Ar / 2 )e-5 A ^V(G)- (4.8) 

Similarly, one can replace ^^(At)! 2 by (/^(O)! 2 in algorithm 2A without affecting its quadratic convergence at very 
small At. We shall refer to this version of the algorithm as 2A0. Algorithm 2B requires two executions of the exact 
algorithm (|3.14|) for similar convergence, which is less efficient. We therefore did not implement algorithm 2B. 

Fig. 1 shows the convergence of algorithm 1A and 1B0 for the chemical potential fi. Both are very linear at small 
At. The calculation is done for 6 = 0.5, £1 = 0.5, and g — 50. This choice corresponds to sizable anisotropy, rotation, 
coupling strength and not close to any particular limit. The calculation uses 64 2 grid points over a 14 2 harmonic 
length square centered on the origin. Changing the grid size to 128 2 only changes the stable results in the fifth or 
sixth decimal place. The ground state wave function is nearly converged by r = 2. The chemical potential shown is 
calculated at t = 10. Note that the linear convergence line for 1B0 fails abruptly at At w 0.15. Since algorithm 2A0 
is just running algorithm 1A first followed by 1B0 at half the time-step size, the convergence failure of 1B0 accounts 
for the failure of algorithm 2A0 near At w 0.3. Both algorithm 1A and 1B0 require an exceedingly small At (< 0.001) 
to produce an accurate value of /i. Even algorithm 2A0 requires At to be « 0.05. 
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V. SELF-CONSISTENT ITERATIONS 



To see the full effect of time-dependent factorization, we must implement IB and 2A in the form (|4.6() and 1|4.7|) 
with self-consistent iterations to determined the GP potential. The required consistency equation is of the form 

4,= J=e-W>, (5.1) 

where b = cArg for some coefficient c, <f> is the unnormalized intermediate wave function prior to the evaluation of the 
potential term, and Z is the constant that normalizes ip. This is needed because we are solving for the GP potential, 
which requires a normalized wave function. Since only the square of the modulus is needed, we solve (|5.1|l as 

x = ^ bx a (5.2) 

where x = \ip(r)\ 2 , a — \<fi(r)\ 2 and 

Z = J e- bx ^a{r) dr » J] e'^a, (Ax) 2 . (5.3) 

i 

It is helpful to view these equations in the discrete forms in which they are actually solved. When necessary, we will 
denote array elements explicitly as Xi — \ip(ri)\ 2 , etc.. 
A naive way of solving (|5.2[) is just to iterate, 

x n+1 = F(x n ) = — l —e- bx ~a. (5.4) 

Z(x n ) 

Starting with xo = 0, the first iteration would produces the normalized a = a/Z, which is a reasonable starting guess. 
Denoting x n = x* + e n , where x* is the exact solution, then 

e n+1 - F\x*)e n + 0{e 2 n ) 

= -bx* (1 - 0(1/N)) E n + ■ ■ ■ 

pa -{cATgx*)e n . (5.5) 

The 0(1 /N) term neglected above is from differentiating 1/Z with respect to X{ (= e~ b ^ Xn ^ i ai/'^2jLi e ^ h ^ Xn ^ a j)> 
where N is the total number of array elements. It is therefore an excellent approximation to regard Z as a constant 
when differentiating with respect to individual array elements. The error propagation H5.5(l explains why this naive 
iteration is undesirable; it will diverge abruptly at large At such that |cArga;*| > 1. 

The self-consistency equation can be solved by better methods, such as Newton's or Halley's iterations. However, 
at large At, the normalized a is not a good enough starting point. Fortunately, l|5.2() has the exact solution 

\w( b 4) (5-6) 



where W(y) is the Lambert W- function 



with series expansion 



25 



b V Z 



We w = y (5.7) 



W(y)^y-y 2 + ^y 3 -ly i + 0(y 5 ). (5.8) 

The series expansion (|5.8|) is not useful for our purpose since its radius of convergence is only 1/e. A better choice is 
the following uniform approximation by WinitzkiSi, 

^)-< 1+S )( 1 -^|l), m 



which has a maximum relative error of 2% (at x f=a 2) in [0,oo]. 
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The normalization constant Z can in principle be determined from 




(5.10) 



but this integral equation is time consuming to solve when W(y) itself has to be computed numerically. A more 
workable scheme is compute compute Z via (|5.3|) from x n , and compute x n+ \ via 1)5. 6|) . The use of l|5.6[) will guarantee 
convergence at all b and At, provided that one can compute W(y) in a simple way. Thus we use the normalized 5, as 
the starting value, compute Z via (|5.3|l . solve for x via (|5.6|) and normalize it to obtain an approximate GP potential 
at all At. We find that this one iteration is sufficient to remove all instability in algorithms 1A and IB; further 
Newton- Raphson iterations produce marginal improvments not worth the additional effort. 

Fig.l shows the convergence of algorithms IB and 2A when the self-consistency condition is approximately satisfied 
by our one W-function iteration. The results are denoted as 1BW and 2AW. The instability in 1B0 and 2A0 no longer 
appears. The convergence of both are linear and quadratic out to large values of At. 

For first and second order algorithms, self-consistent iterations are not needed because At has to be very small in 
order for these algorithms to produce results close to the exact one. If At is small, then one may as well use 1B0 and 
2 AO without wasting time on self-consistent iterations. Self-consistency is a concern only when one is interested in 
enlarging the step-size convergence of higher order algorithms. 



VI. FORWARD FOURTH-ORDER ALGORITHMS 



It is well known from studies of symplectic integrators that factorizations of the form l|4.5|l to (|4.8|) can be generalized 
to higher order in the form 27 i 28 i 29 i 30 i 31 i 32 i 33 i 34 



e -Ar(T+V) = TT p -«,AtT p -biArV 



j-J e -a 4 ArT e - 6i ArV (6 . 1} 

i 

with coefficients {aj, b{\ determined by the required order of accuracy. However, as first proved by ShengS and made 
explicit by Suzuki 3 ^, beyond second order, any factorization of the form (|6.1(l must contain some negative coefficients 
in the set {a^, b{\. Goldman and Kaper— later proved that any factorization of the form Ijtj.lfl must contain at least 
one negative coefficient for both operators. Since a negative time step for the kinetic energy operator will result in an 
unbound and unnormalizable wave function, no such factorization scheme can be used to evolve the imaginary time 
Schrodinger equation, including the Gross-Pitaevskii equation. To go beyond second order, one must use forward 
factorization schemes with only positive factorization coefficientsi2ii 5 rii&. These forward algorithms are currently the 
only fourth order algorithms possible for solving time-irreversible equations with a diffusion kerneUS*ii and have been 
applied successfully in solving the imaginary time Schrodinger equationi2*i£. Omelyand, Mryglod and Folk^iS have 
compiled an extensive list of fourth and higher order symplectic algorithms. However, their sixth- and eight-order 
algorithms contain negative time steps and are not forward algorithms. Recently, one of us have proved^ that while 
sixth order forward algorithms are possible, they require an additional commutator currently not implementable. Thus 
forward algorithms are very unique. Here, we will show that they also yield highly accurate fourth order algorithms 
for solving the Gross-Pitaevskii equation. 

The problem we seek to solve is the ground state of 

H = H X + H y + V(x, y, t) (6.2) 



where V(x, y,r) is the GP potential (|4.4H and 



H x = I p 2 + I(i + ( 5) x 2 -Qxp y , (6.3) 

H y = \pl + \{l-8)y 2 +£lyp x . (6.4) 

The Hamiltonian fundamentally has three operators, which are diagonal in (x,y), (x,p y ) and (p x ,y). If the external 
trapping potential V ex t(x,y) is more general and non- harmonic, we can still write, 

H = H X + H y + V{x, y, t) (6.5) 

but now with 

V(x, y, t) = V ext (x, y)-\{l + 5)x 2 - i(l - 5)y 2 + g\ij>{x, y, t)| 2 . (6.6) 
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The parameter S is then a free parameter associated with algorithm, which we can choose to match the asymmetry 
of Vexti or just set to zero. The crucial points is that, for a rotating trap, harmonic or not, the Hamiltonian has three 
operators diagonal in three separate spaces. By computing the density matrix of 



T = H r 



(6.7) 



exactly via algorithm H3.14|) . we have reduced the Hamiltonian to a two-operator problem. This is a tremendous 
simplification. This simplification is not restricted to harmonic traps, but holds equally for an arbitrary external 
potential. The key point is that the rotating part of the Hamiltonian can be diagonalized regardless of the choice of 
the confining potential. When we diagonalize H x + H y , we generate an inverted harmonic potential in l|6.6|l . which 
must be compensated by the external potential or the GP potential. In the following we will present results only 
for (|6.2[1 . but our algorithm works in the general case of l|6.6|l . We will come back to this point when we discuss 
over-critical rotation in the next section. 

Because implementing T is computational demanding, we must choose a fourth order algorithm with a minimal 
number of T operators. Thus among the many forward algorithms discovered so fani^i^iiLiSiiS, we choose to imple- 
ment only the simplest algorithm, 4A. 



</>(At) 



-iAr\/(Ar) e -lArT g - 



f AtV(At/2) e -|A^ e -lArV(O)^ ) ^ 



with V given by 



V = V 



Ai 



[V, [T, V]} 



Despite the seeming complexity of T as defined by the Hamiltonian (|3.8|) . we have remarkably, 

[V,[T,V]} 

Thus the midpoint effective potential is 



dV\ 2 fdV s2 



dx 



dy 



V(At/2) = <#(Ar/2)|^ + 



, AtV I7<9|V>(At/2)| 2 \2 /«9|?KAt/2)| 2 n2' 



48 



dx 



dy 



(6.8) 



(6.9) 



(6.10) 



(6.11) 



(For the more general case, V is given by l|6.6[) ) The partial derivatives can be computed numerically^ by use of 
finite differences or FFT. Since the FFT derivative converges exponentially with grid size, the use of FFT derivative 
is preferable when the system can be made periodic. In the case with bound state wave functions, this can be done 
by extending the grid size so the the wave function is essentially zero near the grid edge. 

To implement this fourth order algorithm, we first replace ip(Ar/2) and ip(Ar) by ^(O). We will refer to this 
as algorithm 4A00. Its convergence is shown in Fig. 2. We have retained some first and second order results for 
comparison. Aside from its abrupt instability at At w 0.3, its convergence is remarkably flat. All the results at 
At < 0.3 differ only in the fifth decimal place. 

We can also improve the convergence by making the final wave function ?/)(At) consistent with V(At). The results 
for iterating the W-function once as described previously is denoted as 4A0W. By just iterating the W-function once, 
we extended the convergence out to At s» 0.5. Algorithm 4A00 and 4A0W can achieve the result of 2A0 at At nearly 
30 to 40 times times as large. 

To fully implement the time-dependent factorization scheme l|6.8|) and to remove the instability in 4A0W, we evolve 
the midpoint wave function ■0(At/2) from ip(0) by a second order algorithm 2AW and iterate the final wave function 
■(/'(At) for consistency. We denote this algorithm as 4AWW. Its convergence is now smooth, stable and fourth order 
as shown in Fig. 2. 

As pointed out in Ref. Il2t when the eigenfunction converges as 



^ = ^ + O(At"), 

the eigenvalue converges as 

E = E + O(A T 2n ). 
In Fig. 3, we compare the convergence of the GP ground state energy 

E = J r{T+\gm 2 W 2 v/ J H 2 d 2 r. 



(6.12) 
(6.13) 

(6.14) 
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All the fitted lines are of the monomial form 

E - E =CAr n . (6.15) 

Algorithms 1A and 1B0 yielded near-identical quadratic convergence. Algorithm 2AW can be fitted with a fourth order 
monomial as shown. The fit is not perfect because our W-function is only an approximation. Algorithm 4A00 and 
4A0W failed too abruptly to show a smooth trend, but 4AWW can indeed be fitted with an eighth order monomial. 

The computational effort required by each algorithm is essentially that of evaluating the exact algorithm 1|3.14[) . 
which uses 3 2D-FFT. Since 1A, 1B0, and 2A0 all use the exact algorithm once, the second order algorithm 2A0 is 
clearly superior. Algorithms 4A00 requires two evaluations of the exact algorithm plus the gradient potential. The 
gradient potential, if done by FFT, requires 2 2D-FFT. Thus algorithm 4A0n requires 8 2D-FFT, which is 8/3 » 3 
times the effort of algorithm 2A0. Since algorithm 4A00 converges much better than 2A0 at time-steps more than 
three times as large, the class of 4A00 and 4A0W algorithm is clearly more efficient. This efficiency is especially 
evident if higher accuracy is required. The fully implemented algorithms 4AWW use the second order algorithm 
to evaluate midpoint wave function and is therefore w 4 time the effort of 2A0. Looking at Fig. 2, algorithms 
4AWW clearly converge better than 2AW (2A0) even at time-steps four times as large. Note that the first and second 
order algorithms are basically similar, whereas all fourth order algorithm are qualitatively distinct. The second-order 
algorithm is not an order-of-magnitude better than a first-order algorithm, whereas all fourth-order algorithms are 
an order-of-magnitudes better than the second-order algorithm. 

This advantage of fourth order algorithms is cumulative. For example, one can quickly evolve into the ground state 
by use of large time steps. As stated earlier, the GP ground state can be obtained at r = 2. Using algorithm 4A0W 
at At = 0.5, one can get there in four iterations. Algorithm 2A0 would have taken 80 iterations at At w 0.02 . 

To see how these comparisons work out in practice, we give below some timing information. Since running time 
is code and machine dependent, this should be viewed as merely illustrative. The algorithms were programmed in 
Fortran 90 and ran on a 1.2 GH Pentium machine with 0.5 GB of RAM. The IMSL 2D-FFT is used. For the case of 
a 64x64 point mesh, each run of algorithms 1A, 1B0, 2A0, 1BW, and 2AW required 0.0130, 0.0135, 0.0138, 0.0173 
and 0.0203 second, respectively. (Timing is obtained by averaging over 500 run of each algorithm.) Algorithms 1A, 
1B0, 2A0 are indeed comparable, but each self-consistency iteration increases the time by m 0.006 second, which is 
an increase of 40% for 2AW. The time for algorithm 4A00, 4A0W and 4AWW are respectively 0.0394, 0.0459, 0.0630 
second. The time for algorithm 4A00 is approximately three times that of 2A0 and the self-consistency iteration now 
accounts for only an 16% increase. Algorithm 4AWW is slightly more than four times that of 2A0. In the case of 
a 128x128 point mesh, the timing for 2A0, 4A00, 4A0W, 4AWW are respectively 0.0568, 0.1669, 0.1916 and 0.2617 
second. Algorithm 4A00 is accurately three times that of 2A0 and 4AWW is ss 1.5 times that of 4A00. Everything 
scaled up approximately by a factor of four. 

By solving the density matrix of the rotation harmonic oscillator l|6.7|l exactly, we have effected a tremendous 
simplification and has allowed us to derive very compact fourth-order algorithms with excellent large time-step con- 
vergence. They are no more difficult to implement than second order algorithms. If we do not have the exact density 
matrix, then we would have to approximate each occurrence of e _ 2 ATT in ||6.8[) to fourth order, resulting in a much 
more complex algorithm. 

However, by solving the rotating harmonic oscillator exactly, the current algorithms also inherited its limitations. 
As alluded to in Section III, the rotating harmonic trap becomes unstable at O c = y/l — 8. Thus if we are to use 
the exact algorithm Ij3.14|l . we must require f2 < fl c . However, it is known 39-40 that for 8 < 1/5 and g sufficiently 
large, the full GP equation can support "over-critical" rotation in the interval yl — 8 < < \/l + 8. For over-critical 
rotation, one must not prematurely impose the limitation of the rotating harmonic oscillator on the algorithm. 

In light of our previous discussion, we now consider the case of 

V ext (x,y) = 1(1 + A)x 2 + ~(1 - A)y 2 , (6.16) 

and group the Hamiltonian as follow 

H = H x + H y + V(x, y, t) (6.17) 
where H x and H y are defined as before in l|6.3(l . I|6.4[) . and 

V(x, y, t) = ~(A - 8)x 2 - i(A - S)y 2 + g\^(x, y, r)\ 2 . (6.18) 

We thus divorce the deformation parameter S associated with the algorithm, from the physical deformation parameter 
A associated with the trapping potential. If we choose 8 — 0, the algorithm is stable up to £1 = 1, above the physical 
critical value of Cl c = \fl — A. Moreover, since in the Thomas-Fermi approximation the density profile follows the 
shape of the potential, (|6.18|l indeed suggests that the inverted harmonic potential — \Ay 2 can be compensated by 
the GP potential at sufficiently large g, making over-critical rotation possible. 
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VII. CONCLUDING SUMMARY 



In this work we have derived a number of fourth order algorithms and demonstrate their fourth order convergence 
in solving the GP equation in a rotating, anisotropic harmonic trap. These fourth order algorithms, based on forward 
factorization schemes, are the only class of factorization algorithms possible for solving evolution equations with a 
diffusion kernel. Our use of the time-dependent factorization rule provided a systematic way of solving the nonlinear 
GP equations and can be generalized to solve similar nonlinear equation such as the Hartree-Fock and the Kohn- 
Sham equation 4142 . These fourth order algorithms are particularly efficient in solving for the ground state by use 
of large time steps. In constrast to other algorithms, generalizing these algorithms to 3D is very transparent, one 
simply replaces 2D-FFT everywhere by 3D-FFT. Our use of the exact algorithm, which diagonalizes the rotating 
component of the Hamiltonian, is general and can applied to any external trapping potential. This exact algorithm 
also provided insight for understanding over-critical rotation. Physical results obtained by applying these algorithms 
will be presented elsewhere. 
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FIG. 1: Comparing the convergence of first and second-order algorithms in computing the chemical potential of the Gross- 
Pitaevskii equation in a rotating anisotropic trap. The lines are fitted curves to algorithm IA, 1B0 and 2A0 to demonstrate the 
order of convergence of each algorithm. The instability of data points in algorithms 1B0 and 2A0 are removed by the inclusion 
of one self-consistent W-function iteration as indicated by 1BW and 2AW. 
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FIG. 2: The convergence of fourth-order algorithms in computing the chemical potential of the Gross-Pitaevskii equation in 
a rotating trap. Algorithms 4A00 (solid diamonds) and 4A0W (circles) are unstable beyond At m 0.3 and 0.45 respectively. 
Algorithm 4AWW (solid circles) is stable and showed excellent fourth order convergence. 
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FIG. 3: The convergence of various algorithms in computing the ground state energy of the GP equation. Both first order 
results showed near-identical quadratic convergence. The second order result 2AW (asterisks) is fourth order and 4AWW (solid 
circles) is eighth order. Results for 4A0W (circles) cannot be fitted because instability sets in abruptly at At m 0.5. Results 
for 4A00 is similar with instability at ~ 0.3 and is not shown. 



